//## begin module%3C6E5F6D02CD.cm preserve=no
//	  %X% %Q% %Z% %W%
//## end module%3C6E5F6D02CD.cm

//## begin module%3C6E5F6D02CD.cp preserve=no
//	+----------------------------------------------+
//	| SVMBR - A program for training SVMs          |
//	+----------------------------------------------+
//	| * Creator, modeling, coding:                 |
//	|    Marcelo Barros de Almeida                 |
//	|    barros@smar.com.br                        |
//	|    http://litc.cpdee.ufmg.br/~barros/        |
//	| * Coding, improvements, bug fixes:           |
//	|    Bernardo Penna                            |
//	|    bprc@brfree.com.br                        |
//	+----------------------------------------------+
//
//	 Copyright(c) 2002 by Marcelo Barros de Almeida
//	                            All rights reserved
//
//	-[CVS]------------------------------------------
//	$Author$
//	$Date$
//	$Source$
//	$Name$
//	$Revision$
//	------------------------------------------------
//## end module%3C6E5F6D02CD.cp

//## Module: smosolver%3C6E5F6D02CD; Package body
//## Subsystem: <Top Level>
//## Source file: D:\textotese\svmbr\src\smosolver.cpp

//## begin module%3C6E5F6D02CD.additionalIncludes preserve=no
//## end module%3C6E5F6D02CD.additionalIncludes

//## begin module%3C6E5F6D02CD.includes preserve=yes
//## end module%3C6E5F6D02CD.includes

// smosolver
#include "smosolver.h"
//## begin module%3C6E5F6D02CD.declarations preserve=no
//## end module%3C6E5F6D02CD.declarations

//## begin module%3C6E5F6D02CD.additionalDeclarations preserve=yes
//## end module%3C6E5F6D02CD.additionalDeclarations


// Class SMOSolver 










SMOSolver::SMOSolver (REAL cv, REAL tolv, REAL epsilonv, SVM *svmp, Kernel *kernelp, Data *datap, INT chunksizev)
  //## begin SMOSolver::SMOSolver%998688401.hasinit preserve=no
  //## end SMOSolver::SMOSolver%998688401.hasinit
  //## begin SMOSolver::SMOSolver%998688401.initialization preserve=yes
  //## end SMOSolver::SMOSolver%998688401.initialization
{
  //## begin SMOSolver::SMOSolver%998688401.body preserve=yes

	init(cv,tolv,epsilonv,svmp,kernelp,datap,chunksizev);

  //## end SMOSolver::SMOSolver%998688401.body
}


SMOSolver::~SMOSolver()
{
  //## begin SMOSolver::~SMOSolver%.body preserve=yes

	freeSpace();

  //## end SMOSolver::~SMOSolver%.body
}



//## Other Operations (implementation)
INT SMOSolver::run (INT maxi)
{
  //## begin SMOSolver::run%998688365.body preserve=yes

  INT i, j, k;

  // initialize random seed
  SRAND(time(NULL));

  // start to count time
  crono.startCronometer();

  // main loop
  INT numChanged = 0;
  examineAll = TRUE;
  // CHUNK {
  for(i = 0; i < fsSize; i++){
      oldAlphaSubset[i] = alpha[subset[i]];
      fullSet[i] = i;
  }
  oldBiasSubset = bias;
  // }

  cout << "Starting training." << endl << flush;

  // When all data were examined and no changes done
  // the loop reachs its end. Otherwise, loops with
  // all data and likely support vector are alternated
  // until all support vector be found.
  while((numChanged > 0) || (examineAll)) {
    numChanged = 0;
    if(examineAll) {
			heur[0] = 0; 
			heur[1] = 0; 
			heur[2] = 0;	
      cout << "SMO/ALL: " << flush;

	  for(i = RAND(fsSize), j = 0 ; j < fsSize; j++, i++){
				k = i % fsSize;
				numChanged += examineExample(k);
      }
			cout << "IT=" << ++numIter << " CH=" << numChanged << " LS=" << lsSize <<
				     " WS=" << wsSize << " H1=" << heur[0] << " H2=" << heur[1] <<
						 " H3=" << heur[2] << endl << flush;
			//printf("IT=%d CH=%03d LS=%03d WS=%03d H1=%03d H2=%03d H3=%03d\n",
			//			++numIter,numChanged,lsSize,wsSize,heur[0],heur[1],heur[2]);
    } 
		else {
			heur[0] = 0; 
			heur[1] = 0; 
			heur[2] = 0;
			cout << "SMO/KKT: " << flush;
			for(i = 0; i < wsSize; i++){
				numChanged += examineExample(workingSet[i]);
			}
			cout << "IT=" << ++numIter << " CH=" << numChanged << " LS=" << lsSize <<
				     " WS=" << wsSize << " H1=" << heur[0] << " H2=" << heur[1] <<
						 " H3=" << heur[2] << endl << flush;

			//printf("IT=%d CH=%03d LS=%03d WS=%03d H1=%03d H2=%03d H3=%03d\n",
			//				++numIter,numChanged,lsSize,wsSize,heur[0],heur[1],heur[2]);
			if(((wsSize > SUPVEC_PERCENTAGE*fsSize)&&  // chunk
				((numIter % NUM_ITERS_TO_UPDATE_WORKING_SET) == 0)))
					updateWorkingSet();
    }
    if(examineAll == TRUE){
			examineAll = FALSE;
			genWorkingSet(); // generate working set
		}
    else if(numChanged == 0){
			examineAll = TRUE;
			sincSets();// update "lazing set"
		}
  }

	//calcBias();

  // finish time counting
  crono.stopCronometer();
	
  //calcNormW();

  return TRUE;
  
  //## end SMOSolver::run%998688365.body
}

INT SMOSolver::run_chunck (INT maxi)
{
  //## begin SMOSolver::run_chunck%998688366.body preserve=yes

	INT numA = 0, numB = 0;
	INT j, l , k, i;
	PQueue* pqA;
	PQueue* pqB;

	// initing subset: first vectors with 50% of each class
	// data should be randomized
	INT numAMax = data->getnumClassA();
	INT numBMax = data->getnumClassB();
	INT ne      = data->getnumExamples();

	if(numAMax > (chunkSize/2))
		numAMax = chunkSize/2;

	if(numBMax > (chunkSize/2))
		numBMax = chunkSize/2;

	for(i = 0, l = 0, k = 0; i < ne ; i++){
		if((data->gettarget(i) == 1) && (numA < numAMax)){
			subset[k++] = i;
			numA++;
		}
		else 
			if((data->gettarget(i) == -1) && (numB < numBMax)){
				subset[k++] = i;
				numB++;
			}
			else
				supset[l++] = i;
	}


	cache->SetChunk(subset);
	run(maxi);

	// update alpha and bias for remaining set
	
	REAL delta_b = bias - oldBias;
	REAL dalpha;
	
  for(i = 0; i < chunkSize ; i++) {
    k = subset[i];
		dalpha = alpha[k] - oldAlphaSubset[i];
    if (ABS(dalpha) > tol){
      REAL t = data->gettarget(k)*dalpha;

      for(j = 0; j < (ne - chunkSize) ; j++) {
				l = supset[j];
				errorCache[l] += t*cache->evalKernel(k,l);
      }
    }
  }
  
  for(j = 0; j < (ne - chunkSize) ; j++) {
		l = supset[j];
		errorCache[l] -= delta_b;
  }

	// priority queues
	pqA = new PQueue(data->getnumClassA());
	pqB = new PQueue(data->getnumClassB());

	// main loop 
	while(TRUE) {

		// build new subset
		for(i = 0 ; i < ne ; i++){
			if(data->gettarget(i) == 1){
				//if(IS_SV(alpha[i]))
				//	pqA->add(errorCache[i]*10000,i);
				//else
					pqA->add(-data->gettarget(i)*errorCache[i],i);
			}
			else{
				//if(IS_SV(alpha[i]))
				//	pqB->add(errorCache[i]*10000,i);
				//else
					pqB->add(-data->gettarget(i)*errorCache[i],i);
			}
		}

		k = MIN(numAMax,numBMax);

		for(i = 0, j = -1 ; i < k ; i++){
			subset[++j] = pqA->remove(0);
			subset[++j] = pqB->remove(0);
		}

		if(numAMax != numBMax) {
			if(numAMax > numBMax) {
				for(i = 0 ; i < (numAMax - numBMax) ; i++, j++)
					subset[j] = pqA->remove(0);
			}
			else {
				for(i = 0 ; i < (numBMax - numAMax) ; i++, j++)
					subset[j] = pqB->remove(0);
			}
		}

		j = 0;
		while(pqA->getnumElements())
			supset[j++] = pqA->remove(0);

		while(pqB->getnumElements())
			supset[j++] = pqB->remove(0);

		// training
		cache->SetChunk(subset);
		INT oldNumIter = numIter;
		run(maxi);

		// teste convergence
  	if((numIter - oldNumIter) == 1) break;
	
		REAL delta_b = bias - oldBias;
		REAL dalpha;
	
	  for(i = 0; i < chunkSize ; i++) {
			k = subset[i];
			dalpha = alpha[k] - oldAlphaSubset[i];
			if (ABS(dalpha) > tol){
	      REAL t = data->gettarget(k)*dalpha;

	      for(j = 0; j < (ne - chunkSize) ; j++) {
					l = supset[j];
					errorCache[l] += t*cache->evalKernel(k,l);
				}
			}
		}
  
		for(j = 0; j < (ne - chunkSize) ; j++) {
			l = supset[j];
			errorCache[l] -= delta_b;
		}
	
	}
	return true;	
  
  //## end SMOSolver::run_chunck%998688366.body
}

INT SMOSolver::takeStep (UINT i1, UINT i2)
{
  //## begin SMOSolver::takeStep%998688386.body preserve=yes

  REAL y1, alpha1, E1, L, H;
  REAL k11, k12, k22, eta;
  REAL y2, E2, s, a1, a2, alpha2, Lobj, Hobj;
  //REAL f1, f2, L1, H1;
	INT i, j;

  if(i1 == i2)
    return FALSE;

  alpha1 = alpha[i1];
  alpha2 = alpha[i2];
  y1     = data->gettarget(i1);
  y2     = data->gettarget(i2);
  s      = y1*y2;

  E1 = errorCache[i1];
  E2 = errorCache[i2];

  // compute L and H
  if (y1 != y2){
    L = MAX(0, alpha2-alpha1);
    H = MIN(C, C+alpha2-alpha1);
  }
  else {
    L = MAX(0, alpha2+alpha1-C);
    H = MIN(C, alpha2+alpha1);
  }	

  if(L == H)
    return FALSE;
		
  k11 = evalKernel(i1,i1);
  k12 = evalKernel(i1,i2);
  k22 = evalKernel(i2,i2);
  eta = 2*k12-k11-k22;
	
  if(eta <  0 ) {
    a2 = alpha2 - y2*(E1 - E2)/eta;
    if(a2 < L) a2 = L;
    else
      if(a2 > H) a2 = H;
	}
  else {

		
    Lobj = evalObjetive(i1,i2,L);
    Hobj = evalObjetive(i1,i2,H);
		/*
		// eta < 0  : kernel does not obey Mercer's condition
		// eta == 0 : data with bad condition (LD set)
		//Taken from Platt's SMO article (Eqs. 19)
		f1 = y1*(E1+bias) - alpha1*k11   - s*alpha2*k12;
		f2 = y2*(E2+bias) - s*alpha1*k12 - alpha2*k22;
		L1 = alpha1 + s*(alpha2 - L);
		H1 = alpha1 + s*(alpha2 - H);
			Lobj = L1*f1 + L*f2 + L1*L1*k11/2 + L*L*k22/2 + s*L*L1*k12;
		Hobj = H1*f1 + H*f2 + H1*H1*k11/2 + H*H*k22/2 + s*H*H1*k12;
		*/
    if(Lobj > (Hobj+epsilon)) a2 = L;
    else
      if(Lobj < (Hobj-epsilon))	a2 = H;
      else
				a2 = alpha2;
  }

  if(ABS(a2-alpha2) < (epsilon*(a2+alpha2+epsilon)))
    return FALSE;

  // new a1 value
  a1 = alpha1 - s*(a2-alpha2);
  
	/* SVMFU:
  a1 = alpha1 - s*(a2-alpha2);
  if(ABS(a1-alpha1) < (epsilon*(a1+alpha1+1+epsilon)))
    return FALSE;
	*/

  // update bias
  REAL b1, b2, bnew;

	b1 = E1 + y1*(a1 - alpha1)*k11 + y2*(a2 - alpha2)*k12 + bias;
  if ((a1 > epsilon) && (a1 < (C-epsilon)))
    bnew = b1;
  else {
    b2 = E2 + y1*(a1 - alpha1)*k12 + y2*(a2 - alpha2)*k22 + bias;
    if ((a2 > epsilon) && (a2 < (C-epsilon)))
      bnew = b2;
    else
      bnew = (b1 + b2)/2;
  }

  // save bias difference to update error cache and
  // calculate new bias
  REAL delta_b = bnew - bias;
  bias = bnew;
	
  // update error cache
	REAL t1 = y1*(a1-alpha1);
	REAL t2 = y2*(a2-alpha2);


	INT size, *set;
	if(examineAll){
		size = fsSize;
		set = fullSet;
	}
	else {
		size = wsSize;
		set = workingSet;
	}

	REAL * row1 = cache->getKernelRow(i1);
	REAL * row2 = cache->getKernelRow(i2);

  for (i = 0 ; i < size; i++){
    j  = subset[set[i]]; // chunk
    errorCache[j] +=  t1*row1[set[i]] + t2*row2[set[i]] - delta_b;
  }

  alpha[i1]=a1;
  alpha[i2]=a2;
	
 return TRUE;

  //## end SMOSolver::takeStep%998688386.body
}

INT SMOSolver::examineExample (UINT i2)
{
  //## begin SMOSolver::examineExample%998688387.body preserve=yes

  REAL y2, alpha2, E2, E1, r2, tmp, tmax;
  INT i1;
  INT i,j,k;

	INT size, *set;
	if(examineAll){
		size = fsSize;
		set = fullSet;
	}
	else {
		size = wsSize;
		set = workingSet;
	}

	i2 = subset[i2]; // chunk

  y2     = data->gettarget(i2); 
  alpha2 = alpha[i2];
  // calculate distance from (-1),(+1) margems => error
  E2 = errorCache[i2];

	
  // Avoid problems with signals. r2 is negative if point i2
  // is placed between margem (-1),(+1). Otherwise r2 is positive.
  r2     = E2*y2;
	
  // KKT conditions:
  // r2>1  and alpha2==0   (well classified)
  // r2==0 and 0<alpha2<C  (support vectors at margems)
  // r2<0  and alpha2==C   (support vectors between margems)
  // Test the KKT conditions for the current i2 point. 
  //
  // If a point is well classified its alpha must be 0 or if
  // it is out of its margem its alpha must be C. If it is at margim
  // its alpha must be between 0<alpha2<C.
  //
  // If it doens't violate KKT conditions then exit, otherwise continue.
  if(((r2<-tol)&&(alpha2<(C-epsilon)))||
     ((r2>+tol)&&(alpha2>   epsilon))){
    // First the rotine tries to find an i1 lagrange multiplier that
    // maximizes the measure |E1-E2|. As large this value is as bigger 
    // the dual objective function becames.
    // In this first test, only umbounded support vectors will be tested.
    for (i1 = (-1), tmax = 0, i = 0; i < size; i++){
			INT k = subset[set[i]]; // chunk
			E1 = errorCache[k];				
			tmp = ABS(E2 - E1);
			if (tmp > tmax){
				tmax = tmp;
				i1 = k;
			}
    }
    // if a good i1 is found, optimise
    if(i1 >= 0){
			heur[0]++;
      if (takeStep (i1, i2)) return 1;
    }
    // The second heuristic choose any LM that is a SV and tries to
    // optimise
    for(i = RAND(size), j = 0 ; j < size; j++, i++){
			k = i % size;
			i1 = subset[set[k]]; // chunk
			heur[1]++;
			if(takeStep(i1,i2)) return 1;
    }
    // if both heuristc above fail, iterate over all data set 
		if(!examineAll) return 0;
    for(i = RAND(fsSize), j = 0; j < fsSize; j++, i++){
      i1 = subset[i % fsSize]; // chunk
			heur[2]++;
      if(takeStep(i1,i2)) return 1; // TRATAR DIFERENTE !!!!
    }
  }
  return 0;
  
  //## end SMOSolver::examineExample%998688387.body
}

INT SMOSolver::run ()
{
  //## begin SMOSolver::run%998880777.body preserve=yes

	return run(-1);

  //## end SMOSolver::run%998880777.body
}

INT SMOSolver::run_chunck ()
{
  //## begin SMOSolver::run_chunck%998880778.body preserve=yes

		return run_chunck(-1);

  //## end SMOSolver::run_chunck%998880778.body
}

REAL SMOSolver::evalObjetive (UINT i1, UINT i2, REAL lim)
{
  //## begin SMOSolver::evalObjetive%998880788.body preserve=yes

  //chinese algorithm (begin)  
	
  REAL y1 = data->gettarget(i1);
  REAL y2 = data->gettarget(i2);	
  REAL s  = y1*y2;
  REAL ret;
  REAL aa1 = alpha[i1] + s*(alpha[i2] - lim);
  REAL t1  = -y1*aa1/2;
  REAL t2  = -y2*lim/2;

  ret = aa1+lim;
	REAL * row1 = cache->getKernelRow(i1);
	REAL * row2 = cache->getKernelRow(i2);
	
  for (INT i = 0; i < fsSize; i++){
		INT k = subset[i]; // chunk
    if (alpha[k] > 0){
      ret += t1*data->gettarget(k) * row1[i];//evalKernel(i1,i);
      ret += t2*data->gettarget(k) * row2[i];//evalKernel(i2,i);
    }
	}

  return ret;
  
  //chinese algorithm (end)

	//Taken from Platt's SMO article (Eqs. 19)
	/*
  REAL y1 = data->gettarget(i1);
  REAL y2 = data->gettarget(i2);	
  REAL s  = y1*y2;
  REAL a1 = alpha[i1];
  REAL a2 = alpha[i2];
  REAL E1 = errorCache[i1];
  REAL E2 = errorCache[i2];
  REAL f1, f2, LH;
  REAL k11 = evalKernel(i1,i1);
  REAL k12 = evalKernel(i1,i2);
  REAL k22 = evalKernel(i2,i2);
  REAL ret;
  
	// class independent output error without take into account
	// a1 and a2 contributions
	f1 = y1*(E1+bias) - a1*k11   - s*a2*k12;
	f2 = y2*(E2+bias) - s*a1*k12 - a2*k22;
	// value at end of the line segment
	LH = a1 + s*(a2 - lim);
	ret= LH*f1 + lim*f2 + LH*LH*k11/2 + lim*lim*k22/2 + s*lim*LH*k12;

  return ret;
	*/
  //## end SMOSolver::evalObjetive%998880788.body
}

VOID SMOSolver::calcBias ()
{
  //## begin SMOSolver::calcBias%999602348.body preserve=yes
  //## end SMOSolver::calcBias%999602348.body
}

VOID SMOSolver::init (REAL cv, REAL tolv, REAL epsilonv, SVM *svmp, Kernel *kernelp, Data *datap, INT chunksizev)
{
  //## begin SMOSolver::init%1000125177.body preserve=yes

	Solver::init(cv,tolv,epsilonv,svmp,kernelp,datap,chunksizev);

	C       = cv;
	tol     = tolv;
	epsilon = epsilonv;
	svm     = svmp;
	kernel  = kernelp;
	data    = datap;

	bias    = 0;
	oldBias = 0;
	wsSize  = 0;
	lsSize  = 0;
	
	fsSize = chunkSize;/*data->getnumExamples();*/
	try{
		workingSet = new INT  [fsSize];
		lazingSet  = new INT  [fsSize];
		fullSet    = new INT  [fsSize];
		oldAlpha   = new REAL [fsSize];

		for(INT i = 0 ; i < fsSize ; i++){
			workingSet[i] = EMPTY;
			lazingSet[i]  = EMPTY;
			oldAlpha[i]   = 0;
			//fullSet[i]    = i;
		}	
	} 
	catch(...){
		cerr << "Error when allocating memory (SMOSolver::init)." << endl;
		exit(1);
	}

  //## end SMOSolver::init%1000125177.body
}

VOID SMOSolver::freeSpace ()
{
  //## begin SMOSolver::freeSpace%1000125178.body preserve=yes

	Solver::freeSpace();

	delete oldAlpha;
	delete workingSet;
	delete lazingSet;
	delete fullSet;

  //## end SMOSolver::freeSpace%1000125178.body
}

VOID SMOSolver::sincSets ()
{
  //## begin SMOSolver::sincSets%1000125179.body preserve=yes

	INT i, j, k, l, m;
	REAL delta_b = bias - oldBias;
	
  for(i = 0; i < wsSize ; i++) {
    k = workingSet[i];
		m = subset[k];
    if (ABS(alpha[m] - oldAlpha[k])>tol){//alpha[k] != oldAlpha[k]) {
      REAL t = data->gettarget(m)*(alpha[m]-oldAlpha[k]);

			REAL * row = cache->getKernelRow(m);

      for(j = 0; j < lsSize ; j++) {
				l = lazingSet[j];
				m = subset[l];
				errorCache[m] += t*row[l];
      }
    }
  }
  
  for(j = 0; j < lsSize ; j++) {
		l = lazingSet[j];
		m = subset[l];
		errorCache[m] -= delta_b;
  }

  //## end SMOSolver::sincSets%1000125179.body
}

VOID SMOSolver::genWorkingSet ()
{
  //## begin SMOSolver::genWorkingSet%1000125180.body preserve=yes

	wsSize = 0;
	lsSize = 0;
	oldBias = bias;

  for(INT i = 0; i < fsSize ; i++){
		INT k = subset[i]; // chunk
		oldAlpha[i] = alpha[k]; // chunk
		if(IS_SV(alpha[k]))//((alpha[i]>epsilon)&&(alpha[i]<(C-epsilon))) // chunk
			workingSet[wsSize++] = i;
		else
			lazingSet[lsSize++] = i; 
  }
  //## end SMOSolver::genWorkingSet%1000125180.body
}

VOID SMOSolver::updateWorkingSet ()
{
  //## begin SMOSolver::updateWorkingSet%1000125181.body preserve=yes
	cout << "Updating working set..." << endl << flush;
	sincSets();
	genWorkingSet();
  //## end SMOSolver::updateWorkingSet%1000125181.body
}

// Additional Declarations
  //## begin SMOSolver%3B87FD8602AC.declarations preserve=yes
  //## end SMOSolver%3B87FD8602AC.declarations

//## begin module%3C6E5F6D02CD.epilog preserve=yes
//## end module%3C6E5F6D02CD.epilog
